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Abstract 

Generalized-ensemble algorithm and diffusion theory have been combined in or- 
der to compute the dynamical properties monitored by nuclear magnetic resonance 
experiments from efficient and reliable evaluation of statistical averages. Replica- 
exchange Monte Carlo simulations have been performed with a C-peptide analogue 
of ribonuclease A, and Smoluchowski diffusion equations have been applied. A fairly 
good agreement between the calculated and measured ^H-NOESY NMR cross peaks 
has been obtained. The combination of these advanced and continuously improv- 
ing statistical tools allows the calculation of a wide variety of dynamical properties 
routinely obtained by experiments. 
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1 Introduction 



In the study of protein folding, a crucial step is the understanding of secondary structure 
formations. The construction of a-helices and /5-sheets from disordered structures and 
their interconversion is mainly driven by hydrophobic effects combined with dispersive 
interactions and intramolecular hydrogen-bond formations ^1 12] . Salt bridges and strong 
electrostatic interactions can either compete with this driving force or assist the secondary 
structure formation, by introducing long-range average distance constraints. The latter 
can be influenced by the characteristics of the solvent and/or ions eventually enclosed in 
pre-folded configurations. 

In order to address these complications related to the sequence variability in real pro- 
teins' secondary structures, many experiments are performed on relatively short peptide 
fragments that in vitro adopt conformations similar to those revealed in the entire pro- 
teins. These fragments, therefore, may be considered as independent modules jHj of the 
original protein These experiments allow a great simplification, replacing most of the 
protein matrix with the solvent, while keeping in the molecule the relevant interactions. 
On the other hand, short peptides visit many conformations, and thus the usage of sta- 
tistical tools and computer simulations is required for the interpretation of experimental 
data. 

Computer simulations in canonical ensemble based on Boltzmann weight factor, how- 
ever, tend to get trapped in states of energy local minima, and it is very difficult to obtain 
accurate statistical averages even for small peptide systems. Generalized-ensemble algo- 
rithms are based on artificial, non-Boltzmann weight factors and perform random walks 
in potential energy space, which efficiently alleviates the multiple-minima problem (for 
a recent review, see Ref. After a single production run, one can calculate accurate 
canonical-ensemble averages for a wide range of temperatures jSj. Because of the very con- 
struction of the generalized-ensemble algorithms, however, the information of molecular 
dynamics is lost, and only static average values can be obtained by these methods. 

Nuclear magnetic resonance (NMR) has become an attractive technique because it al- 
lows one to monitor both structural statistics (average distances) and molecular dynamics 
(stochastic rotation of vectors) at an atomistic level and in different solvent /temperature 
conditions. These features are, however, intimately linked together in the experimental 
data, particularly in those data routinely measured to obtain structures, i.e., the ^H- 
NOESY cross peaks (CPs). 

In order to calculate dynamical properties from NMR experiments, a diffusive model 
has been designed and casted in a Smoluchowski diffusion equation. This equation has 
been solved by matrix expansion methods and by using a mode-coupling approximation [01 
[7]. This procedure allows the description of the time-correlation functions (TCFs) that 
govern NMR experiments through the computation of a suitable set of configurational 
averages. 

In this Letter we propose to combine the above two powerful computational ap- 
proaches: generalized-ensemble algorithm and diffusion theory, which allows rigorous 
calculations of molecular dynamics inferred by NMR experiments. 

A C-peptide analogue of ribonuclease A is here considered. The C-peptide of ribonucle- 
ase A is one of the smallest peptides that is known to form a-helix conformations and has 
been extensively studied by circular dichroism (CD) [Hj and NMR jHl CHj spectroscopies. 
The peptide has also been studied by computer simulations J^^C^]- Nevertheless, a 
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quantification of the conformational population and of its effect on the experimental data 
was not possible and the standard methods used to analize NMR relaxation data cannot 
be applied ^3] . We employed one of the commonly used generalized-ensemble algorithms, 
replica-exchange Monte Carlo (REMC) ^Hj, to calculate statistical averages and various 
dynamical quantities were successfully obtained by the diffusion theory. 

In section 2 the two mehotds are summarized and computational details are given. In 
section 3 the results of the statistics and the NMR CPs are discussed in a unique frame. 
In section 4 conclusions and perspectives are presented. 



2 Methods 

2.1 Replica- exchange method 

We first briefly review the replica-exchange method (REM) JHI (see, for instance, Refs. f?! 
E] for details). 

The system for REM consists of M non-interacting copies (or, replicas) of the original 
system in the canonical ensemble at M different temperatures Tm {m = 1, ■ ■ ■ , M). We 
arrange the replicas so that there is always exactly one replica at each temperature. 
Then there is a one-to-one correspondence between replicas and temperatures. Let X = 
| - ■ ■ , xj^l, ■ ■ - j stand for a "state" in this generalized ensemble. Here, xj^ stands for the 
state of the i-th replica (at temperature Tm)] the superscript i and the subscript m in 
xj^l label the replica and the temperature, respectively. Each state x]^ is specified by the 
coordinates (and momenta pW) of all the atoms in replica i. 

A simulation of REM is then realized by alternately performing the following two 
steps. Step 1: Each replica in canonical ensemble of the fixed temperature is simulated 
simultaneously and independently for a certain MC or MD steps. Step 2: A pair of 
replicas, say i and j, which are at neighboring temperatures and T„, respectively, are 
exchanged: X = |- ■ ■ , ■ ■ ■ , x]^\ ■ ■ - j — > X' = {■ ■ ■ , a^m, " " " > a;[^V ' '}• The transition 
probability of this replica exchange is given by the Metropolis criterion: 

w{X^X')^w{x^^ |a;|^])=min(l,exp(-A)) , (1) 

where 

A = (/5„-/5„)(E(g[^])-E(g['l)) . (2) 

Here, E (^g'''^ and E (^q^^^^ are the potential energy of the i-th replica and the j-th replica, 
respectively. In the present work we employ Monte Carlo algorithm for Step 1. When the 
potential energy depends on temperature as in the present case (see Eq. (fT^ below), we 
should use the following A instead of that in Eq. Q [IE] (see also [TTI]): 

A = /3m [E (g[^]; Tm) - E (g^; Tm)) - i^n (e [q^^; T„) - E (g^; T„)) . (3) 

A random walk in "temperature space" is realized for each replica, which in turn 
induces a random walk in potential energy space. This alleviates the problem of getting 
trapped in states of energy local minima. 

The canonical expectation value of a physical quantity A at temperature Tm (fn = 
1, • • ■ , M) can be calculated by the usual arithmetic mean as follows: 

< ^ >T.= E ^ (^l^(^)) ' (4) 

sim i—\ 
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where Nsim is the total number of measurements made at temperature Tm- Note that 
the above summation is taken over different rephcas i {i = 1,- ■ ■ , M) that happens to 
correspond to the fixed temperature at the moment of measurement. The expectation 
values at any intermediate temperature can, in principle, be calculated by the multiple- 
histogram reweighting techniques ^Zj, but in this article we limit our discussions to the 
above M temperature values. 



2.2 Calculation of NMR parameters 

In this subsection we briefly review recent advances of diffusion theory [01 Ej that are 
applied to the calculation of the ^H-NOESY NMR cross peak intensities (CPs). 

The CP intensities Oh,k at mixing time between spin h and spin k can be calculated 
by the following equation: 

R 

Oh,k (tm) = exp (-a tm) sinh (A tm) , (5) 

where Rh,k, ^, and a are all functions of the following spectral densities (see Refs. [HI EI 
for details): 

oo 

Jh,k{uj) = 2 1 COS (ut) TCFh^k (t) dt . (6) 



Here, TCF^ k is a time-correlation function of 2nd-rank tensor components of the vectors 
r/i joining proton h and proton k. These TCFs at temperature T have the form 



TCF{t)= Y: {[D^Zi^it)) /r{t)'][D^SAm) /m']). (7) 

M=-2 

where i^^o irreducible spherical tensors [21], and f2 and r are the direction and the 
modulus, respectively, of the given H-H vector involved in the NOESY CP (subscripts h 
and k are henceforth suppressed for clarity). 

In order to separate the effect of the modulus from the orientation of the unit vector 
(direction) in the above TCP, we also calculate the orientational TCP from 

TCFoit) = E {[D^m))] [DZim)])^ = P2(cos[e(t)]) , (8) 

M=-2 

where P2 is the 2nd-order Legendre polynomial. The orientational mobility can be de- 
scribed by the correlation time r that is the integral of TCFq'- 



T 



J TCFo (t) dt . (9) 



The mode-coupling diffusion (MCD) theory of the dynamics of a biological macro- 
molecule in solution is adopted for the computation of the above TCPs of Eqs. ((Tj) and 
(jH}. The MCD approach [221 12S] can be briefly summarized as follows. Given a poly- 
mer of Na beads of friction coefficients Q and coordinates rj, connected by Nf, bonds 
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(Ij, i = 1, Nf,), the dynamics of each variable U, is regulated by the operator L, adjoint 
to the diffusion Smoluchowski operator D: 

- = LI; L=Y: [V.D,, V, - (V^^U/kBT) D,, V,] , (10) 

where U is the potential energy of the beads as a function of the bead coordinates, fc^ is 
the Boltzmann constant, T is the absolute temperature, and I is the 3 x A^;, dimensional 
array containing all the bond vectors Ij. 

By expanding the conditional probability (solution to the Smoluchowski equation) in 
a complete set of eigenfunctions of L, the time autocorrelation function (TCF) of any 
coordinate-dependent dynamic variable with zero average f{t) may be expressed in the 
standard form 

(/W/(0))t = Y.{f^^)T{^^f)Texp{-\t), (11) 
i 

where — Aj and ipi are respectively the eigenvalues and the normalized eigenfunctions of 
the operator L: 

L^i = -\iiPi . (12) 

This eigenvalue equation becomes a matrix equation, with the matrix elements being 
equilibrium averages at temperature T (see Refs. jHl E] for details). In the present work, 
we use the replica-exchange Monte Carlo method for the calculation of these averages (see 
Eq. (H). 



2.3 Computational details 

The configurational statistics of a C-peptide analogue jH] with the amino-acid sequence 
AETAAAKFLRAHA and uncharged N- and C- termini have been simulated by replica- 
exchange Monte Carlo method. Residue His 12 was protonated in order to better match 
the NMR experimental conditions of pH 5.2. Other charged residues were Glu 2~, Lys 
T"*", and Arg lO"*". The number of atoms in the model was 195. 

The total "potential energy" function E{q] T) that we used is the sum of the confor- 
mational energy term of the solute Ep{q) and the solvation free energy term Esol{(1]T) 
for the interaction of the peptide with the surrounding solvent: 

E{q;T) = Ep{q) + EsoL{q;T) . (13) 

The parameters in the conformational energy as well as the molecular geometry were taken 
form ECEPP/2 The sigmoidal, distance-dependent dielectric function of Ref. P3] 

was used. 

The solvation free energy that we used is given by a linear combination of the solvent- 
accessible surface area (SASA) Ai of each non-hydrogen atom i: 

EsoLiTo)=J2a,A, , (14) 

i 

where cxj are the proportionality constants, Tq = 298 K, and the dependence on the 
coordinates q is now suppressed. The temperature dependence of the solvation free energy 
was taken into account, following the prescription in Ref. 



T / T 

Esol{T) = —Esol{To) + Hsol{Tq) I 1 - — ) - Csol(^o) 



Tln(^) + Tn-T 



Tr 



(15) 
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where Hsol and Csol are enthalpy and heat capacity, respectively [22] • 

The SASA was calculated by the computer code NSOL j^Zj. The computer code 
KONF90 was used, and MC simulations based on REM were performed. In each 
Monte Carlo sweep all the independent dihedral angles except for the peptide-bond di- 
hedral angles cu, which were fixed at 180 degrees, were updated once and the Metropolis 
test was performed for each update. The number of changeable torsion angles was then 
55 and the maximum torsional change was ±180 degrees. 

For REMC we used 10 replicas. The corresponding temperatures were 200, 233, 
276, 317, 370, 432, 504, 588, 686, and 800 K. These temperatures were chosen to span 
the temperature range between 200 and 800 K and to contain the temperature of the 
NMR experiments (276.15 K). The initial conformations were randomly generated. The 
replica exchange was tried every 20 MC sweeps. For an optimal performance of REMC 
simulations, the acceptance ratios of replica exchange should be sufficiently uniform and 
large (say, > 10 % ). The acceptance ratio was indeed found to be in the range 14-24 %, 
and we observed that each replica underwent an unbiased random walk in the potential 
energy space (and that each temperature underwent a random walk in the replica space). 
After 1,100,000 MC sweeps of equilibration, the REMC production run of 3,000,000 MC 
sweeps for each replica was made. The configurations were stored every 10 MC sweeps 
for data analyses. This amounts to 300,000 configurations for each temperature (or each 
replica) . 

As for the diffusion equation, the first step is to approximate atoms or group of atoms 
in the molecule as friction points. In the present model, friction points were located on 
55 heavy atoms (beads) among the total of 195 atoms. The friction was computed by 
using Stokes' law with stick boundary conditions with Stokes' radii obtained summing 
the accessible surface area to a spherical probe of radius (ASAO) of the atoms grouped 
in the bead |2H1- The Stokes' radii ranged from 0.09 nm (Cq, in all residues that include 
Hq only) to 0.24 nm (the last portion of Arg 10 side chain, that includes N^, C^, N^i, N^2 
and all the bonded hydrogen atoms). The water viscosity 7] was 0.001 Pa s and it was 
assumed independent of temperature. 

The second step is to assess the convergence of the used basis set in solving the 
eigenvalue equation for the adjoint of the diffusion operator. Both the long-time sorting 
procedure (LTSP) [SHI and the maximum correlation approximation (RM2-II basis set of 
MCA) j7] were applied to select the most important terms of the infinite mode-coupling 
basis set. The results coupling the five Ist-order lowest-rate modes for Ist-rank variables 
up to the 2nd-order for 2nd-rank variables (MCA with e = 5 in the notations of Ref. [7j) 
were almost identical with the 2nd-order LTSP using up to 600 basis functions (data not 
shown). Therefore, the MCA basis set built with e = 5 (240 basis functions) was used for 
all of the following calculations. 

No significant changes in the ten lowest-rate Ist-order /Ist-rank relaxation modes were 
observed reducing the number of configurations from 300,000 down to 10,000. Therefore, 
for the calculations of the statistical averages required to solve the diffusion equation and 
to compute the TCFs at each temperature, 10,000 of the 300,000 recorded configurations 
were used. 
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3 Results and Discussion 



The C-peptide is relatively rich in hydrophobic residues (Ala, Leu, and Phe) and, there- 
fore, in a water environment is expected to be mainly in a-helix conformation. On the 
other hand, the presence of charged residues, namely, Glu 2~ (near the N-terminus), 
Lys 7"*", Arg 10^, and His 12~^ (near the C-terminus), will have significant effects on the 
conformational states of the peptide. In Fig. ^ the average a-helicity as a function of 
residue number (or, probability of each residue being in the a-helix state) is shown for 
four different temperatures (276 K, 370 K, 504 K, and 686 K). Here, we considered that a 
residue is in the a-helix state when the backbone dihedral angles (0, ip) fall in the range 
(—70 ± 30°, —37 ± 30°), and Eq. (0)) was used to calculate the average helicity. At the 
lowest temperature (T = 276 K) among the four, residues 4-12 are in helical state (espe- 
cially, residues 6-11 are completely helical), which is exactly the same location of a-helix 
as found in the corresponding structure from the X-ray experiments of the entire ribonu- 
clease A jHI]]. Helicity decreases as the temperature is raised because of the increased 
thermal fluctuations. High helicity persists up to T = 504 K (especially in residues 6-11), 
and finally at the highest temperature among the four (T = 686 K), extended helical 
conformation ceases to exist. Note also that the N-terminus is not helical even at 276 K 
in agreement with the NMR experiments jH]. 

Fig. 121 shows the behaviour of the average total helicity (or, average total number of 
helical residues) as a function of temperature. As was observed in Fig. 1, the residues 
near the N-terminus are rarely in helical state, and this is the reason why even at the 
lowest temperature (T = 200 K) the average total number of helical residues is only 
about 8. The slight decrease of helicity between 370 K and 500 K is due to residues 6, 7, 
and 12 that lose a-helical population, while the further decrease beyond 500 K involves 
the demolition of the remaining a-helix in residues 8-11. The disorder in configurational 
statistics beyond 600 K still keeps part of the electrostatic interactions characterizing the 
low-energy structures (see discussion below), thus representing a stiff disordered polymer 
segment. Note that our statistics show that high total helicity (> 50 %) persists as high as 
T ^ 500 K, while the experiments observe high helicity only near T = 273 K [Hlini- This 
shift in helix-coil transition temperature is presumably due to the fact that our energy 
functions including the solvent model are not accurate enough to reproduce the absolute 
temperature dependence of experiments. As discussed below in detail, our simulation 
results around T = 400 K best reproduce the NMR experiments (which were conducted 
at T = 276 K). 

In Fig. inithe end-to-end distance distribution is shown at the same four temperatures 
as in Fig. ^ Here, the end-to-end distance is defined to be the distance between N of 
Ala-1 and O of Ala- 13. At the two lower temperatures (T = 276 K and 370 K) we observe 
three peaks in the distributions, which suggests that there exist three groups of similar 
conformations. We refer to the three groups as Groups 1, 2, and 3 from left to right in the 
Figure. Representative conformations from each group (the lowest-energy conformation 
in each group) are also shown in Fig. IHl All these three groups of conformations have 
a common a-helix structure in residues 5-11. The end-to-end distance is about 1.5 nm, 
2 nm, and 2.5 nm for Group 1, Group 2, and Group 3, respectively. Conformations in 
Group 1 are characterized by a salt bridge between side chains of Glu 2" and Lys and 
a bend towards the N-terminus so that the end-to-end distance is the shortest among the 
three groups. Conformations in Group 2 are charecterized by two salt bridges between 
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Glu 2~ and Lys 1^ and between Glu 2~ and Arg 10"^. Note that this group has the most 
similarity to the X-ray structure jHO]! which also has the Glu 2- Arg 10 salt bridge. The 
backbone root-mean-square distance of the Group 2 structure in Fig. 3 from that of X-ray 
experiments is 0.12 nm, while those of Group 1 and Group 3 are 0.32 nm and 0.19 nm, 
respectively. Finally, conformations of Group 3 also have a salt bridge between Glu 2~ 
and Lys 7^ . However, the N-terminus is pushed away from the a- helix and the structure 
is rather extended. 

As is shown in Fig. 3, the highest populated peak at the lowest temperature (276 
K) corresponds to Group 1. As temperature increases (370 K), this peak decreases in 
population and the peak of Group 2 increases, thus showing an increased stability of a 
longer a-helical segment (characterized by a larger end-to-end distance). At the third 
temperature (504 K) only a single peak (of Group 2) exists, which suggests that the 
"native-like" structure (i.e.. Group 2) is the most stable among the three groups. At 
the highest temperature (686 K) we have a single peak at a different end-to-end distance 
(about 2.2 nm), which corresponds to a coil structure. 

The above results can be explained in terms of microscopic interactions. Glu 2 is 
involved in salt bridges with Lys 7 and, less frequently, with Arg 10. These salt bridges 
bend the N-terminus toward the short a-helical region in residues 6-12 and tend to make 
the whole molecule more compact. These interactions are more efficient at the low tem- 
peratures (276 K). Increasing the temperature, the observed salt bridges become less 
populated and the bending is more frequenly released, thus allowing an increase of the 
molecular extension and, eventually, a more extended a-helical region including the Ala- 
Ala-Ala sequence. These extended a-helical configurations are actually present in the 
statistics and are characterized by an end-to-end distance of about 2.5 nm. Beyond 500 
K, the hydrophobic interactions mainly responsible for the a-helical stability are no more 
effective. On the other hand, the salt bridge between Glu 2 and Lys 7 is still contained 
in the statistics (although it is weaker). 

We now study the effects of this conformational distribution and of its temperature 
dependence on the NOESY CPs. The NOESY experiment was performed at the tem- 
perature of 276.15 K, pH 5.2, in the static magnetic field corresponding to vi^H) = 500 
MHz, and using a mixing time of 400 ms. These experimental CPs are reported in refer- 
ence jH]. Unfortunately, the experimental results are reported in arbitrary units and do 
not allow a quantitative estimate of global dynamical effects, but it is possible to analyze 
the behaviour of CPs for different proton pairs in the molecule. 

The experimental CPs are reported for three types of proton pairs: HN(i)-HN(i+l) 
(referred to as NN CPs, hereafter), Ho(i)-HN(i-|-l) (AN CPs, hereafter), and some long- 
range CPs involving side chains (LR CPs, hereafter). First, the NN CPs were computed 
at four different temperatures that were used in the REMC simulation (i.e., 276, 370, 504 
and 686 K). The above analysis of several statistical quantities, such as the helicity for 
each residue at each temperature (Fig. H]) showed that the molecular model undergoes 
a helix-coil structural transition with a transition temperature estimated between 450 K 
and 600 K (Fig. El). The real sample undergoes the same transition at lower temperatures; 
the experimental NOESY CPs are not observable even at room temperature. In Figs. Ufa) 
and llfd) the behaviour of the NN and AN CPs with temperature is shown. In both NN 
and AN data sets, negative CPs are present for the lowest temperature. The experimental 
NOESY pulse- sequence is phase-sensitive and CPs of opposite sign with respect to the 
diagonal peaks are not observed at the experimental conditions. Explanation of the change 
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in sign of the CPs is not trivial. The explanation in terms of a global r, which is the 
inverse of the unique relaxation rate usually assumed to govern the 2nd-rank rotation 
of the given H-H distance vector, is here not meaningful, because of the high flexibility 
that allows many rates to play a significant role in each spectral density. The change in 
sign of CPs can be related to the change in the whole 2nd-rank rate spectrum that is 
obtained by the calculations (data not shown); at the lower temperature the spectrum is 
characterized by several gaps that are progressively smoothed by increasing temperature. 
It is expected that for the highest temperature the spectrum be almost a smooth function 
of the relaxation mode as in a polymer random coil. The intermediate situations, where 
the internal kinetics is faster, but the molecule is still characterized by groups of internal 
modes separated by rate gaps, can produce the change of sign in CPs and CPs relatively 
high in magnitude, as in the experiments. On the other hand, the large CPs calculated 
at high temperatures (e.g., the results at T = 686 K in Fig. ^ can be related to the 
limitations of the basis set construction in the MCA approximation [i.e., e=5). 

Even if the orientational mobility of the H-H unit vector and the H-H average distances 
are coupled in the TCF of Eq. the behaviour of CPs with temperature can be partially 
understood in terms of the correlation time (Eq. Q) and of the average H-H distances. In 
Figs, m^b) andlU^e) the correlation times for the NN and AN unit vectors are respectively 
plotted, and in Figs.|3fc) and|3ff) the average moduli of the same proton pairs are plotted. 
Correlation times globally decrease as is expected. The decrease in orientational rigidity 
of the low-temperature helical region occurs up to the highest temperature where the 
difference in orientational rigidity between residues 1-3 and 4-11 is not significant, as it 
is expected for a non-structured molecule. On the other hand, it must be noticed that 
the distances behave differently: The differences between distances in the two regions 
disappear at the highest temperature, but the distances are still small enough to give 
contribution to the CPs, especially the AN CPs. Therefore, the structural information 
contained in CPs must be searched in the behaviour of both sets of CPs, and a separate 
analysis of AN and NN CPs may be misleading. 

The orientational mobility is more sensitive to temperature than the average distances, 
the latter being more stable. However, it is evident that the approximation of assuming 
the same orientational behaviour for all the H-H vectors along the sequence does not hold; 
at 370 K an increase of r from 50 ps to 300 ps moving the H-H vector from the N-terminus 
deep into the a-helix (Phe 8) can be observed. 

As a consequence of the analysis of the temperature behaviour of CPs, the statistics 
between 370 K and 504 K can be considered to reproduce the experimental conditions 
at T = 276.15 K. A shift in the order-disorder transition temperature with respect to 
experiments is expected and always occurs in other simulated systems like liquid crystals. 

In Fig. El CPs at four selected temperatures are compared with the available experi- 
mental data. All calculated and measured CPs are scaled by a unique factor in order to 
have the AN CP of Ala 3, which is the largest experimental CP available, equal to 1. The 
increase of NN CPs from the N-terminus to the helical region (residues 6-12) is reproduced 
by calculations at T = 370 K, with the largest deviations beyond Phe 8 (Fig. Efa)). The 
significant decrease of AN CPs beyond residue 3 is qualitatively reproduced (Fig. Efb)); 
CPs beyond residue 5 are slightly smaller than those found in the experiments and the 
largest deviations from experiments is found for Ala 4. This latter deviation is also found 
in the ^J(HN-Ha) coupling constants (data not shown). The Coupling constants have 
been computed using the Karplus equation = 1.9 — 1.4 cos(0') + 6.4 cos(0') where 
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0' is the HN-N-Cq-Hq, dihedral angle. Residue Ala 4 presents the largest deviations from 
experiments, thus implying that the torsional state in the region 4-5 is not well captured 
by the simulation. 

It must be noticed that NN CPs are better reproduced at 370 K, while AN CPs seem 
to be better reproduced at 432 K. This suggests that a qualitative reproduction of both 
sets of CPs could be achieved at an intermediate temperature. 

The computed LR CPs are all weak; the strongest one (H^(Thr 3)-HN(Ala 5) at 
T = 370 K) is 1.2 x 10"'^ compared to the experiment where it is found to be about the 
same order of magnitude of the NN CPs. These low CPs are caused by the large average 
distances, because the orientational correlation times of LR H-H unit vectors are in the 
range 100-200 ps, therefore only slightly smaller than those of NN H-H unit vectors in 
the helix (about 200 ps, see Fig. life)). On the other hand, the H^(3)-HN(5) average 
distance is 0.5 ±0.1 nm, which is too large to produce even a weak CP. However, the Glu 
2-Lys 7 salt bridge that was suggested by these long range NMR constraints is found very 
stable up to T = 432 K in the model (0.4 nm), while it becomes 0.5 nm at the highest 
analyzed temperature. Therefore, even if the NMR structural details are not reproduced 
in terms of H-H distances, the global features that are responsible for the most populated 
conformations are contained in the model up to the temperature where the comparison 
between experimental and computed NN and AN CPs is qualitatively good. 

4 Conclusions 

In this work two recent advances in statistical mechanics have been combined together in 
the study of the statistics and dynamics of a small peptide, the C-peptide of ribonucle- 
ase A. The replica-exchange Monte Carlo simulation has been used to sample molecular 
configurations in the canonical ensemble at several temperatures in the range 200-800 
K. This method, together with other generalized-ensemble algorithms, has the advan- 
tage of allowing the overtaking of energy barriers, connecting low temperature and high 
temperature trajectories. Diffusion theory in the form of the Smoluchowski equation for 
the conditional probability governing the stochastic time evolution of intramolecular seg- 
ments' orientation, has then been used to model the orientational correlation functions 
and to compute the ^H-NOESY NMR cross peaks that are experimentally available. 

The combination of generalized-ensemble statistics and diffusion theory, frequently 
updated by technical progresses making both methods more robust and efficient, allows 
the direct calculation of NMR data and other dynamical properties, thus closing the gap 
between theoretical or computational models and experiments. 
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Figure Captions 

• Figure ^ Average helicity as a function of residue number at four temperatures: 
T = 276 K (solid line), T = 370 K (dashed line), T = 504 K (dotted line) and 
T = 686 K (dotted-dashed line). 

• Figure 1^ Average total helicity as a function of temperature. 

• Figure IHl Distribution of the end-to-end distance at four temperatures: T = 276 
K (solid line), T = 370 K (dashed line), T = 504 K (dotted line) and T = 686 
K (dotted-dashed line); the arrows identify the peaks in terms of representative 
conformations. Besides backbone, the side chains of Glu 2, Lys 7, and Arg 10 are 
also shown. The conformations were drawn with MolMol |31j . 

• FigurelU NN (a-c) and AN (d-f) CPs, correlation times r, and average H-H distances 
as functions of residue number: T = 276 K (squares), T = 370 K (circles), T = 504 
K (triangles) and T = 686 K (diamonds). 

• FigureEl NN (a) and AN (b) CPs as functions of residue number: experimental data 
(filled squares) and calculated data at T = 276 K (squares), T = 370 K (circles), 
T = 432 K (triangles), and T = 504 K (diamonds). 




Figure 1: La Penna et al. 
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Figure 2: La Penna et al. 
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Figure 3: La Penna et al. 
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Figure 4: La Penna et al. 
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Figure 5: La Penna et al. 



